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^ . Abstract 

Cn ' We investigate effects of hardness of primordial binaries on whole evolution of star clusters by means 

J>^| of A^-body simulations. Using newly developed code, GORILLA, we simulated eleven A^ = 16384 clusters 

C^ . with primordial binaries whose binding energies are equal in each cluster in range of 1 — SOOfcTo, where 

l.SfcTo is average stellar kinetic energy at the initial time. We found that, in both soft (< SfcTp) and hard 

(> SOOfcTo) limits, clusters experience deep core collapse. In the intermediate hardness (10 — lOOfcTo), the 

CO ■ core collapses halt halfway due to an energy releases of the primordial binaries. The core radii at the halt 

can be explained by their energy budget. 
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1. Introduction 

Recently, observational informations concerning on the binary systems in the globular cluster have been accumulated 
from photometric observations of the eclipses (e.g. Mateo 1996), spectroscopic observations (e.g. Albrow et al. 2001), 
low-mass X-ray binaries (see Liu et al. 2007), and also the color-magnitude diagrams (e.g. Rubenstein et al. 1997). 
Davis et al. (2008) also constrain the binary fraction in NGC 6397 by the method of the color-magnitude diagrams, 
and review the binary fractions in galactic globular clusters. 

The binaries in the globular cluster play important roles for its dynamical evolution, since the binaries work as 
energy source through interactions with other stars or binaries (Heggie 1975). Even if the globular cluster has no 
binary at initial, the binaries are formed in due course through three-body encounter, and energy generated from 
these binaries halts core collapse of the globular cluster (Henon 1975). 

Furthermore, if the globular cluster contains a nonnegligible fraction of primordial binaries, the dynamical evolution 
could be affected. Goodman and Hut (1989) (hereafter GH89) first showed theoretically difference of evolutions with 
and without primordial binary. GH89 estimated that the core of the cluster with primordial binaries at the halt of core 
contraction is larger than that of the cluster without primordial binaries by order of magnitudes. GH89's estimate is 
based on the model that the cluster core stop contracting when energy generated by the primordial binaries is balanced 
with energy outflowing from the inner region of the cluster to the outer region through two-body relaxation. McMillan 
et al. (1990; 1991) first performed iV-body simulations of clusters with primordial binaries. They clearly showed that 
the cluster cores stop contracting at larger cores than those without primordial binaries. By means of Fokker-Planck 
model, Gao et al. (1991) investigated post-collapse evolution of cluster with the primordial binary. They showed that 
the cores of the clusters continues to contract slowly after rapid core contraction and that the gravothermal oscillations 
occur after the several ten half-mass relaxation time of the slow core contraction. 

In order to set a population of the primordial binaries in the cluster, several parameters concerning the binary and 
its distributions, such as mass fraction, and distribution of binding energies and eccentricities, need to be specified. 
Although these parameters should be derived from theories of star and cluster formations or be limited observationally, 
sufficient informations have not yet been provided at present. Only a few studies investigated the effect of these 
parameters on the evolution. Heggie and Aarseth (1992) performed A^-body simulations of clusters with two different 
mass fractions, 6 and 12 %, of primordial binaries, and Heggie et al. (2006) investigated the evolutions of clusters 
with 0— 100 % mass fraction of primordial binaries by means of A^-body simulations. Vesperini and Chernoff (1994) 
(hereafter VC94) extended the model of GH89, taking into account the core mass fraction of primordial binaries and 
the distribution of the binding energies of the binaries. 

In this paper, we focus on the dependence on binding energies of the primordial binaries. By means of A^-body 
simulations, we systematically investigate the dynamical evolution of clusters with different initial hardness of binding 
energies. We set the distribution of binding energies of primordial binaries by delta function, (5(a; — i?bin.o)i where 
i?bin,o is the initial binding energies of the primordial binaries and i?bin,o = 1,3, 10,30, 100, SOOfcTo. Here, 1.5fcTo is 
the average kinetic energy of stars in the cluster at the initial time. In previous simulations, the distribution of the 
binding energy of the primordial binaries are usually fixed at uniform distributions in log_Ebin.o (see McMillan et al. 
1990; Gao et al. 1991; Heggie and Aarseth 1992; Heggie et al. 2006). 

We found that the evolutions of the cores are different according to the distribution of the binding energies. If 
primordial binaries consist of softer binaries, the binaries are disrupted before core collapse through binary-single and 
binary-binary encounters, and do not affect the core evolution. If primordial binaries consist of harder binaries, the 
binaries escape from the cluster. This is because energy generated by the harder binaries per one encounter is larger 
than the cluster potential. They also do not affect the core evolution and the cluster exhibits deep core collapse. In 
intermediate range, the binaries efficiently heat the clusters, and the core collapses of the clusters halt halfway. These 
behaviors are consistent with the theoretical estimate by VC94. 

These behaviors were not shown by previous works (McMillan et al. 1990; Gao et al. 1991; Heggie and Aarseth 
1992; Heggie et al. 2006; Fregeau, Rasio 2007; Trenti et al. 2007) who followed the dynamical evolution of clusters with 
primordial binaries. This is mainly because the previous works fixed the distributions of the binding energies of the 
primordial binaries at uniform distributions in logi^bin.Oj fmd their fixed distributions include binaries with hardness 
in soft, intermediate, and hard ranges. 

The structure of this paper is as follows. We describe simulation methods in section 2. In section 3, we present 
results of simulations and investigate the effect of the hardness of primordial binaries on the dynamical evolution of 
clusters. In section 4, we summarize this paper. 

2. Simulation methods 

In this section, we describe simulation methods. In section 2.1, we show initial conditions of cluster models. We 
perform A^-body simulations of the clusters by means of GORILLA, which is a newly developed A-body simulation 
code for star clusters (Tanikawa, Fukushige 2009, hereafter TF09). In section 2.2, we outline GORILLA. 
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2.1. Initial conditions 

By means of A^-body simulations, we simulate the dynamical evolution of 13 clusters with point-mass particles and 
without an external tidal field. The cluster models are shown in table 1. Each of 11 clusters have equal- mass stars 
and contains primordial binaries with equal binding energy. Among these cluster models, binding energies of the 
primordial binaries, -Ebin.Oj a-nd mass fractions of the primordial binaries, /b,o- are different, as shown in the second 
and third columns of table 1, respectively. The number of the primordial binaries iVb,o is also shown in the fourth 
column of table 1. The remaining two clusters are reference models as soft and hard limits; an A'' = 16384 equal-mass 
cluster model without primordial binaries, and a cluster model in which all binaries are replaced by stars with double 
mass in the /b,o = 0.1 models. The fifth and sixth columns in table 1 are the mass fraction of the double mass stars, 
/d,Oi and the number of the double mass stars, iVd,o, respectively. 

In all models, we use Plummer's model to generate the initial distribution of both single stars and center of mass 
of primordial binaries in the clusters. The eccentricity distribution of the primordial binaries is thermal distribution, 
f{e)de = 2ede. The other orbital elements of the primordial binaries, such as the inclination, the longitude of the 
ascending node, and the argument of pericenter with respect to the clusters, and the phase are distributed at random. 

We adopt iV-body standard units (Heggie, Mathieu 1986), such that G = M = —E = 1, where G is the gravitational 
constant, M is the total mass of the cluster, E is the total energy of the cluster not including the internal binding 
energy of the primordial binaries. For clusters with 16384 stars, IfcTo = 1-0 x 10~^. 

2.2. GORILLA: an N-hody simulation code for globular clusters 

In order to perform A^-body simulations of the above clusters, we use an TV-body simulation code for the star cluster, 
GORILLA (TF09). In GORILLA, the orbits of cluster stars are integrated with a fourth-order Hermite scheme with 
individual timestep (Makino, Aarseth 1992). The timesteps of stars are quantized by power of 2 (McMillan 1986) in 
order to be adjusted to GRAPE, a special-purpose computer designed to accelerate A^-body simulations, where in our 
simulations we use GRAPE-6/6A (Makino et al. 2003; Fukushige et al. 2005). Additionally, the relative motions of 
two stars relatively isolated from other stars are approximated as Kepler motions. 

In order to determine timestep of each star in the fourth-order Hermite scheme with individual timestep, we use 
the following criterion. 



lal a(2) -I- a(^) ^ 



^t=iV-^ 9: (1) 

" |a(i)||a(3)| + |a(2)|' 

where a is the acceleration of each star, a'"' is n-th order derivative of a, and rj is an accuracy parameter. We set the 
accuracy parameter ?y = 0.01. We also use the following timestep criterion for startup, 

At^^^j^, (2) 

1^ I 

where 77s is a startup accuracy parameter. We set the startup accuracy parameter r]s = 0.0025. 

In GORILLA, the relative motions of binaries are approximated as Kepler motion if the binary components (stars 
k and I) satisfy either of two conditions as follows. 

Isolation conditions (A) 

1. |r3-rcm,fe/| > arapo,fc/ 
Isolation conditions (B) 

1. efe/>0.95 

Here, r^ is the position of the nearest star originated from the center of mass, rcm,ki, of stars k and I, Tapo.fc; is the 
separation between stars k and I at the apocenter, rroi,fe/ is the separation between stars k and I at a given time, and 
Cki is eccentricity. 

Condition (A) expresses whether the third star is separated enough compared to apoccntric distance of a binary. 
Condition (B) expresses whether the third star is separated enough compared to the instantaneous distance of the 
binary when the binary is highly eccentric. The above conditions are illustrated in figure 1. Using dimensionless 
quantities a (called apocentric parameter) and /3 (called pericentric parameter), we give criterions of the isolation. 
We show the apocentric and pericentric parameters adopted in each cluster model in table 2, as well as the accuracy 
parameters, 77 and r]s. 

We set both the isolation conditions (A) and (B) in the following reason. Consider a highly eccentric binary which 
does not satisfy isolation conditions (A). Unless rapo.fcz of the binary arc changed, the binary is not regarded as an 
isolated binary even if its pericentric distance is very close. At the pericenter, the accuracy of the orbital calculation 
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is drastically decreased due to rioi./t;. Owing to isolation conditions (B), the binary is regarded as an isolated binary 
around its pericenter. 

We apply the isolation not only to binaries, but also to unbound two stars and hierarchical triple systems. Such 
unbound two stars are sufficiently isolated from other stars and close to each other. If the two stars satisfy the second 
condition of isolation conditions (B), they are regarded as isolated stars. If both timesteps of the two stars are less 
than 2~^^ in standard A^-body units, they are regarded as close stars. In the hierarchical triple system, a binary and 
one star orbit around each other. The binary in the hierarchical triple system is called inner binary. If the inner 
binary is approximated as point mass, we may regard the inner binary and the other star as a binary. Such binary is 
called outer binary. If the hierarchical triple system is in isolation, both the inner and outer binary satisfy isolation 
conditions (A). 

3. Results 

3.1. Accuracy of N -body simulations 

Figure 2 shows energy errors as a function of simulation time in all cluster models. They are all within 1 % (~ 0.0025) 
of the total energy not including the internal binding energy of the primordial binaries at the initial time. 

It seems that the energy errors 1 % are relatively large for studies of the evolution of cluster cores whose energies 
are 1 — 10 % of the total energies of the clusters, excluding the total binding energies of the binaries. However, the 
energy errors do not much affect the core evolutions. We found that large part of the energy errors attributes to the 
binding energy of the binaries in the clusters when we use GORILLA (TF09). 

We also found that the energy errors do not affect properties of the binaries, such as distribution of binding energies 
of the binaries. The energy errors are at most 0.0025 as seen in figure 2. On the other hand, the total binding energies 
of the binaries increase by about 0.2 for models No-binary and Double, and all /b.o models (see figure 4). The energy 
errors are only ^--^ 1 % of the increased amount of the total binding energy in each model. 

3.2. Core evolution and binary properties 

We first see the core radii, re, and the half-mass radii, r^, in the /b,o = 0.1 models. Figure 3 shows the time evolution 
of the core radii and half-mass radii of six /b,o = 0.1 cluster models, and models No-binary and Double. We calculate 
the core radii as in Casertano and Hut (1985) with the modifications described in McMillan et al. (1990). We calculate 
the core and half-mass radii at each time unit, and average these radii over 10 time units. 

In models No-binary, IfcTg — 0.1, and SfcTg — 0.1, the clusters experience deep core collapse, and gravothermal 
oscillations occur. The core radii at the halts of the core collapse are 0.002 — 0.004. In models lOfcTg — 0.1, SOfcTg — 0.1, 
and lOOfcTo — 0.1, the core collapse stops halfway, and the cores contract more slowly. The core radii at the halts of 
the core collapse are 0.05 — 0.1. In models SOOfcTo — 0.1 and Double, the clusters also experience deep core collapse. 
In model Double, gravothermal oscillations occur. The core radii at the halts of the core collapse are 0.005 — 0.02. 
Among the clusters that experience deep core collapse, the times when the core collapse stops and core bounce occurs 
are different. In models No-binary and IfcTo - 0.1, i ~ 3400, in model BfcTo - 0.1, i -- 4700, and in models SOOfcTo - 0.1 
and Double, t ~ 1700 - 2100. 

Next, we see the evolution of binary properties throughout the rest of this subsection. Figure 4 shows in the thick 
curves the increase of the total binding energy of the binaries, AiJbin,tot(i), in the /b,o = 0.1 models and models 
No- binary and Double. The arrows indicate the times when the core collapse stops. The increase AE'bin tot(i) is given 
by 

Afb(t) Wb(0) 

ASbi„,tot(i) = Y. ^bi„.^(i)- Yl ^bi„,.(0), (3) 

i i 

where -Ebin,j(0 i^ the binding energy of i-th binary at time i, Nb(t) is the number of the binaries at time t including 
binary escapers, and Ai^^bin^tot (i) corresponds to energy released by all the binaries. We can see that, in models No- 
binary, IfcTo — 0.1, and Double, the binaries do not release energy until the core collapse stops. In model 3kTo — 0.1, 
the binaries release energy from t ^ 1000. In the other models, the binaries release energy from t ~ 0. After the 
core collapses in models No-binary, IfcTo — 0.1, SfcTo — 0.1, and Double, three-body binaries release energy (discussed 
below). 

Figure 4 shows in the thin curves the time evolution of the total kinetic energy of escapers, i?osc,tot(i), in the 
/b,o = 0-1 models and models No-binary and Double. The escapers are defined as stars, regardless of single stars, 
binaries, or hierarchical triple systems, satisfying both conditions as follows. 

(a) The sum of the kinetic and potential energy of the single star (or the center of mass of the binary or hierarchical 

triple system) is positive. 

(b) The distance between the star and the center of the cluster is more than 40 length units. 
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We can see that, in models No-binary, l/cTp — 0.1, and Double, the total energy of the escapers is small just before the 
halts of the core collapse, such that £'csc,tot(0 ^ 1 x 10~^. In model SOOfcTo — 0.1, the thin curve is almost overlapped 
with the thick curve, which indicates that the escapers carry away almost all energy released by the binaries. In the 
other models, i?csc.tot(0 ^^ the halts of core collapse is larger than those of models No-binary, IfcTo — 0.1, and Double 
by an order of magnitude. 

Figure 5 shows the time evolution of the number of binaries, TVb) in the /b,o = 0.1 models and models No-binary 
and Double. For model Double, the number of the double mass stars is also plotted. In each panel, the thick curve 
shows the number of binaries (or double mass stars) within the cluster, and the thin curve shows the total number of 
binaries (or double mass stars) including escapers. The arrows indicate the times when the core collapse stops. 

We can see that, in model No-binary, the binaries increase after deep core collapse. These binaries are the three-body 
binaries. In models IfcTg — 0.1 and 3fcTo — 0.1, the numbers of the binaries rapidly decrease before deep core collapse. 
After deep core collapse, the total numbers of the binaries including binary escapers increase. In these models, the 
three-body binaries are also formed. In models lOfcTo — 0.1, SOfcTo — 0.1, lOOfcTg — 0.1, and 300/cTo — 0.1, the numbers 
of the binaries monotonically decrease. In models No-binary, IfcTo — 0.1, and 3fcTo — 0.1, the three-body binaries are 
formed after deep core collapse. When the simulations are finished, the numbers of the three-body binaries are 29, 23, 
and 15, and the numbers of the escapers of the three-body binaries are 22, 19, and 15 in models No-binary, IfcTo — 0.1, 
and 3fcTo — 0.1, respectively. In model Double, the three-body binaries are also formed after deep core collapse. The 
total number of these binaries including the escapers is 27 at t = 5000. The number of the binaries composed of the 
two double mass stars is 25. Two binaries arc composed of one double and one single mass stars. Nearly all binaries 
are composed of the double mass stars. The number of the binaries within the cluster is 2, both of which are composed 
of the two double mass stars. In this model, the number of the double mass stars monotonically decreases. 

Figure 6 shows the number of binaries (N^) in each logarithmic bin of binding energies of the binaries (-Ebin) for 
the /b,o = 0-1 rnodels at the time indicated in each panel. All the binaries in the clusters at those time are counted, 
and are primordial binaries. There arc no three-body binaries. In all models, the peaks around the initial binding 
energies, E'bin, can be seen. The distributions of the binding energy spread towards the larger sides. 

In figure 7, the binding energies are shown as a distance from of the cluster center for all the binaries in models 
IfcTo -0.1, 3fcTo-0.1, lOfcTo-0.1, 30fcTo-0.1, lOOfcTo-0.1, and 300fcTo-0.1 at the time indicated in the panels. Ah 
the binaries in the clusters at those time are counted, and are primordial binaries. The dashed lines show the half-mass 
radii and twice the core radii at the time. Within the half-mass radii, the distributions of the binding energies are 
greatly changed from those at the initial time in all the models except model 300fcTo — 0.1. The distributions of the 
binding energies in these models are similar. The distributions center on about lOOfcTo, and range from more than 
lOfcTo to 300fcTo. On the other hand, outside the half-mass radii of these models, the distributions of the binding 
energies arc little changed from the initial time. In model 300fcTo — 0.1, even inside the half-mass radii, the distribution 
of the binding energies is not changed so much. 

In model IfcTo — 0.1, there are little binaries between 2rc and r^. The binaries do not increase between 2rc and rh 
despite of mass segregation, unlike the core. Furthermore, the binaries are disrupted through binary-single encounters 
even between 2rc and rh. The large semi-major axes of the IfcTo binaries result in frequent binary-single encounters 
despite of stellar density in this region smaller than that in the core, and the softness of the binaries results in the 
destruction of the binaries through the encounters despite of the average stellar kinetic energy in this region smaller 
than that in the core. 

Figure 8 shows the time evolution of the mass fraction, /b, of the binaries inside the core radii (upper curves) and 
half- mass radii (lower curves) in the /b.o = 0.1 models and model No-binary. For model Double, mass fraction, /d, of 
the double mass stars are shown. In model IfcTo — 0.1, both the fractions inside the core and half- mass radii decrease 
from 0.1 to 0.02 until the deep core collapse occurs. In model 3fcTo — 0.1, the fraction inside the core radius increases 
up to 0.3 at i ^ 1000, and decreases down to ~ 0.04 at the time when the deep core collapse occurs. In models 
lOfcTo — 0.1, SOfcTo — 0.1, lOOfcTo — 0.1, 300fcTo — 0.1, and Double, the fraction inside the core radii increases until the 
core collapse stops. After that time, the fractions stop increasing. 

Figure 9 shows the time evolution of the mean kinetic energies of the single stars, i?kin,avc,s, inside the core radii 
(solid curves) and the half-mass radii (dashed curves) in the /b,o = 0.1 models and models No-binary and Double. In 
models No-binary, IfcTo — 0.1, 300fcTo — 0.1, and Double, at the moment of core collapse, the mean kinetic energies of 
the single stars in the cores increase several times more than those at the initial time. On the other hand, they nearly 
keep constant during the evolution in models 3fcTo — 0.1, lOfcTo — 0.1, SOfcTo — 0.1, and lOOfcTo — 0.1. 

Figure 10 shows the time evolution of the mean kinetic energy of the binaries, -Ekin.ave.b, inside the core radii in the 
./b.o = 0.1 models. We do not show those inside the half- mass radii, since they are nearly the same as inside the core 
radii. The fiuctuation is large in model IfcTo — 0.1, since the number of the binaries is small. Except model Ifc'To — 0.1, 
the mean kinetic energy of the binaries in the core keep nearly constant during the evolution. In model IfcTo — 0.1, 
the mean kinetic energy largely increases at the deep core collapse. 



6 A.Tanikawa and T.Fukushige PASJ [Vol. , 

3.3. Interpretation 

On both softer (models No-binary, IfcTo — 0.1, and 3kTo — Q.l) and harder (models SOOfcTb — 0.1 and Double) 
hardness, the clusters undergo deep core collapse. On the other hand, in the intermediate hardness (models lOfcTg — 0.1, 
SOfcTo — 0.1, and lOOfcTo — 0.1), the clusters exhibit shallower core collapse. The depth of the core collapse depends 
on the amount of energy heating core generated by the primordial binaries. The larger the amount of energy is, the 
shallower core collapse becomes, and vice versa. 

The amount of the energy heating the core depends on whether the primordial binaries become harder or not through 
binary-single and binary-binary encounters, and whether the single stars and binaries heated by such encounters are 
ejected or not from the clusters. Whether the binaries become harder or not depends on whether the binding energy of 
the binaries is larger or not than a critical energy i?crit,H, which, in this case, corresponds to the average kinetic energy 
of the surrounding single stars, £'kin,avo,s (Heggie's law: Heggie 1975). If E'bin ^ EaitM "^ -£^kin,avo,s, the binaries are on 
average destroyed through binary-single encounters. Therefore, the binaries cannot heat the core. If -Ebin ^ -E'crit.Hi 
the binaries become harder and harder through series of binary-single encounters, and then heat the core. Since most 
binary-single encounters occur in the core, wc should adopt the average kinetic energy of the single stars in the core. 
In Plummer's model, which is the initial condition of our models, the average kinetic energy of the single stars in the 
core is about 2.5fcTo, which can be seen in figure 9. Therefore, 

^crit3 ^ -Ekin.avc.s ^ 2.5fcT, (4) 

where l.SfcT is the average stellar kinetic energy in the whole cluster at a given time. 

Whether the single stars and binaries heated through encounters are ejected or not depends on whether the kinetic 
energies transformed from the binding energies of binaries are larger or not than the potential depth of the whole 
cluster. The increase of the binding energy at single binary-single encounter, A_Ebin, is on average AE'bin — 0.4i?bin 
(Heggie 1975). According to conservation of momentum, two thirds of the energy released by the binaries go to the 
single star, and the rest goes to the binary, on average. The condition to eject both the single star and binary is 
Aii^bin/3 > 2m|$c|, where m is the mass of the single mass stars, and <I>c is the potential of the core. If we define 
a critical energy for ejection as i^^crit.E, 0.4iJcrit.E/3 > 2m|$c|- In Plummer's model, to|$c| — lOfcT. Therefore, the 
ejection occurs when 

-Bbin>^crit,E-^™|*c| = 150fcr. (5) 

In summary, whether primordial binaries can heat the core or not are different among three ranges of hardness divided 
by two critical hardness Ecrit.n and £^crit,E; (a) £^bin ^ £^crit,H' (b) £^crit,H ^ £'bin ^ £^crit,E, and (c) E'bin ^ ^'crit.E, which 
are illustrated in figure 11. Here and hereafter, binaries with smaller binding energy than E'crit.H are called soft, and 
those with larger binding energy than i?crit,H are called hard. If -Ebin ^ ^crit,H (soft range, (a)), the primordial binaries 
are destroyed through encounters, and can not heat the core. If i?crit,H ^ E\^in ^ -Bcrit.E (intermediate hard range, (b)), 
the primordial binaries become harder and harder and continually heat the core. If E'bin ^ £'crit,E (super hard range, 
(c)), the primordial binaries release their binding energies to the surrounding single stars and binaries, but the single 
stars and binaries are ejected from the cluster, and then the primordial binaries can not so much heat the core. In 
the following sub-subsections, we see the simulation results for each range of hardness. 

3.3.1. Soft range 

In the soft range, -Ebin ^ £'crit,H('~^ 2.5fcT), the primordial binaries are destroyed through encounters, and cannot 
heat the core. We can see, in figure 5 (the second left panel), the number of binaries for model IfcTo^ 0.1 rapidly 
decreases, and, in figure 4 (the second left panel), primordial binaries do not release energy from the beginning. After 
then, their evolutions are almost identical to the case without primordial binaries, and lead to deep core collapse. 
Figure 3 and 4 (the first and second left panels) show that the core evolution and energy generations are very similar 
between models No-binary and IfcTg — 0.1. 

Model 3/i;To — 0.1 exhibits a mixed behavior between soft and intermediate hard ranges. While the number of the 
binaries rapidly decreases, shown in figure 5 (the third left panel), the binaries continually release energy, shown in 
figure 4 (the third left panel). A population whose binding energy is harder than initial value can be seen in figure 6 
(the second left upper panel). Although that heating is not so large to stop core collapse, it makes the time to the 
core collapse longer, compared to the case without primordial binaries, shown in figure 3 (the first, second, and third 
left panels). Note that at the very beginning the mass fraction of binary in the core increase, shown in figure 8 (the 
third left panel), due to mass segregation of single star and binary, and sinking of binary to the core. After t = 1000, 
the mass fraction of binary in the core turns to decreasing, since the binaries are destroyed in the core. 

3.3.2. Intermediate hard range 

In the intermediate range, £'crit,H('~-" 2.5fcT) ^ -Ebin ^ £'crit,E(~ 150fcT), the primordial binaries continually heat the 
core, and becomes harder and harder. We can see in figure 4 (the fourth, fifth, and sixth left panels) the primordial 
binaries release the binding energy, and in figure 6 (panels lOfcTo^ 0.1, SOfcTo — 0.1, and lOOfcTo — 0.1), populations 
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whose binding energy become harder from initial values can be seen. The continually released heat halts the core 
collapse halfway. 

3.3.3. Super hard range 

In the super hard range, -Ebin ^ -E'crit,E(^ ISOfcT), the kinetic energy transformed from binding energy of primordial 
binaries through encounters is so large to be ejected from whole cluster immediately, and the primordial binaries cannot 
heat the core. Figure 4 (the second right panels) shows the released binding energy in total (in the thick curves), and 
those of escapers (in the thin curve) are very close in model SOO/cTq — 0.1, which means almost all released energy 
from the primordial binaries is brought away from the cluster by the escapers. Therefore, the primordial binaries 
can neither heat the core, nor stop core collapse. Their evolutions become similar to the case in which binaries are 
replaced by the double mass stars, which is shown in figure 3 (first and second right panels). 

3.3.4. Theoretical estimate 

We theoretically estimate core size at the halt of core contraction, assuming that the energy outflowing from the 
inner region of the cluster to the outer region at each unit time, dEh/dt, is balanced with the energy provided for the 
core through binary interactions at each unit time, dEc/dt, such as 

dEh dEc ,„, 

The argument here and hereafter is based on VC94. 

The energy outflowing from the inner region of the cluster to the outer region at each unit time, dE^/dt, is given by 

dEi^^\E\_^(UGM^ 

dt 7ti.h 7 ^rhrh ' 

where G is the gravitational constant, M is the total mass of the cluster, and 7 is a numerical coefficient relating the 
energy outflow rate to the half-mass relaxation time. The half-mass relaxation time i,h is expressed as 

_ 0.138Mi/2rh3/2 
^*~Gi/2^1og(0.47V)' ^' 

where and rh is the average mass of its stars (Spitzer 1987). We use \E\ ^ 0.2GM'^/rh, which comes from virial 
theorem. \E\ ^ GM'^ /Ar^, where r^ is the virial radius, and rh ~ O.Sr^ in Plummer's model. 

The energy provided for the core by binary interactions at each unit time, dEc/dt, is expressed as 

dEc ,^, G^TO^ 9 / . , N 

-T^ = XVc nc" Abs + ^bb 9 

dt Ws.c 

where Vc is the core volume, Vg^c is one dimensional velocity dispersion of single stars in the core, n^ is the total number 
density of the single stars and binaries in the core, and Abs and ^dbb are, respectively, dimensionless efficiency factors 
for energy provided for a cluster through binary-single and binary-binary encounters. The core volume, Vc, is given 
by K; = (47r/3)rj?. We introduce a free parameter A, which is nearly equal to unity, and actually we adopt A = 0.45 
as discussed later. The free parameter A is the fraction of the core volume where most energy is released through 
binary-single and binary-binary encounters in the unit of the core volume. 

Substituting equations (7) and (9) into equation (6), we express the ratio of the core to half- mass radii at the halt 
of core collapse as 

where /b,c is mass fraction of binaries in the core, Vh is the one dimensional half-mass velocity dispersion, and 
\E\ ~ 3Mwh^/2 i.e. 3wh^/2 ^ GM/bri^. Here, we define the core radius as 



9«c2 _ /9z;,,c2(2-/b,c)2 



AnGpc y IdTrGmric 

^(2 — /b,c)/2 and pc ~ 27nric/(2 — /b,c) are, respectively, the average velocity dispersion and mass density 
in the core. 

The dimensionless efficiency factors of binary-single and binary- binary encounters are indicated by Abs and Abb- 
We derive these factors below, although the derivations are nearly the same way as VC94. This is because two points 
are different from theirs. One is the unit of energy. In this paper, fcT indicates the average stellar kinetic energy in 
the whole cluster. In VC94, however, kT indicates the average stellar kinetic energy in the cluster core. The other is 
the numerical factor of the cross section of binary-binary encounters, 5, which is first seen in equation (19). 

The dimensionless efficiency factors are, respectively, expressed as 



Abs 
and 



2(1 -/b,c) 



2 - /b., 



/b.. 



2 - /b., 
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fix)[g{x)hix)+g'ix)]dx, (12) 



^bb = 77 7^ 'j— 1 / f{xi)f{x2)G{xi,X2)H{xi,X2)dxidx2, (13) 



2 V 2 — /bx 

where x, Xi, and 0:2 are the bmding energy of the binaries in the unit of kT. The function f{x) is the distribution 
function of the binding energies of the binaries in the core. The functions g{x), g'{x) and G{xi,X2) are, respectively, 
the dimensionless hardening rates of the binary with the binding energy x which is not destroyed in a sea of single 
stars, the binary with the binding energy x which is destroyed in a sea of single stars, and the binary with the binding 
energy xi in a sea of binaries with the binding energy X2- The functions h{x) and H{xi,X2) are, respectively, the 
efficiency ratios of heating of the core to hardening of the binaries at each interaction between a single star and a 
binary with the binding energy x, and that between binaries with the binding energies xi, and X2- The h{x) and 
H{xi,X2) become less than unity when single stars and binaries involved with encounters are ejected from the cluster 
immediately after the encounters. When the binaries are destroyed through binary-single encounters, h[x) = 1, since 
they are hardly ejected. 

The dimensionless hardening rate g{x) is expressed as 

ff(x) = 1.66(g)'Yi?b,(.T/C,A)AdA, (14) 

where i?bs is the dimensionless rate of the interactions that the binary with the binding energy x hardens to the 
binding energy (1 + A)a:: in a sea of single stars, and C(= 1.7) is a correction factor in order to set the unit of the 
binding energy to be kTc{= 1.7kT), which is one dimensional kinetic energy of single stars in the core. The integral of 
the dimensionless rate i?bs over A is described in equation (49) of Heggie and Hut (1993), and obtained as a function 
of the binding energy of the binary in the unit of kTc- 
The dimensionless hardening rate g'{x) is expressed as 



X 



-1 



g'{x) = -lM[-) RU^/C), (15) 

where R[^^ is the dimensionless rate of the interactions that the binary with the binding energy x is destroyed in a 
sea of single stars. The dimensionless rate i?(jg is described in equation (5.12) of Hut and Bahcall (1983), and also 
obtained as a function of the binding energy of the binary in the unit of kTc- 
The dimensionless hardening rate G{xi,X2) is expressed as 

G{xi,X2) = 1m(^^^^] R,,b{xuX2)A, (16) 

where i?bb is the dimensionless rate of the interaction that the binary with the binding energies xi hardens to the 
binding energy (1 + A)(2;i -I-X2) in a sea of binaries with the binding energy X2 when xi > X2, and the binaries with 
the binding energy X2 are destroyed. We consider only binary-binary interaction which results in the destruction of 
the softer binary, since Mikkola (1983a; 1983b; 1984a; 1984b) showed that the binary-binary interaction not involving 
the destruction of one binary has small contribution to the heating of the cluster. The dimensionless rate is averaged 
over A, since the number of binary-binary scattering experiments is much smaller than that of binary-single scattering 
experiments. The average of A in binary-binary interactions is about 0.5. 
The dimensionless rate i?bb is expressed as 

Rhh{xi,X2) = — 5 / vaj{v)dv, (17) 



i^s,c 



where a is semi- major axis of the binary with the binding energy a;i -I- 2:2 , v is the relative velocity between binaries, j (v) 
is the distribution of the relative velocity between the binaries, and a is cross section of the binary-binary interactions. 
Assuming that the velocity distribution of the single stars and binaries in the core is isotropic Maxwellian, and 
equipartition is achieved: fb,c = (l/v2)^s,c! where Wb,c is one dimensional velocity dispersion of the binaries, the 
relative velocity dispersion between the binaries v^^'b.c is Ws.c, and j{v) is expressed as 



s.c 



JM=(z) — e^P(-^^)- (18) 



The cross section of the binary-binary interaction a is expressed 

'y^S—^ : (19) 
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where i5bin,2 is the binding energy of the softer binary, expressed as X2mVg ^/C, and S* is a dimensionless coefficient. 
The cross section a is derived, based on equation (2.7) in Gao et al. (1991). The dimensionless coefficient S depends 
on the relation between the binding energies of the binaries, xi and X2, as follows: 



5 = 



25.2 
15.9 



{Xi '-X2) 

{xi >a;2) 



(Gao et al. 1991). 

In summary, the dimensionless rate i?bb is expressed as 



C 



X2 



and the dimensionless hardening rate G{xi,X2) is expressed as 



G(a;i,X2)= 0.3985 



Xi +X2 
X2 



(20) 



(21) 



(22) 



We describe the forms of h{x) and H{xi,X2)- When we obtain h{x), we assume that at every binary-single encounter 
the binary increases its binding energy, E'bin, by 0.4iJbin, and two thirds of the increment go to the kinetic energy of 
the single star and the rest to the kinetic energy of the center of mass of the binary. The increment is the average 
value over all the binary-single encounters, which has been obtained by Heggie (1975). The single star will be ejected 
when E\y/m\^c\ > 15/4. The binary will be ejected when E-t,/m\^c\ > 15. Since the ejection results in the mass loss of 
the core, the binding energy of the core decreases, i.e. the core is heated. The amount of the heating is m-|$c| when 
a single star is ejected, and 277i|<I>c| when a binary is ejected. Therefore, we can express h{x) as 

15TO|(f>rl 



h{x) 



m|(f>r 



Ai^bs 
3to|$c| 

AE, 



if a; < 

ff ^ <x<lh- 

4 kT 

li x>15 



kT 



(23) 



bs 



kT 



We set m|$c| = lOfcT as described in section 3.3. Then, h{x) is expressed as 
if a; < 38 ; 
ff 38 < .T < 150; 

iix> 150. 



h{x) 



1, 
1 



25 



(24) 



For H(xi,X2), we assume that at every binary-binary interaction, 0.5(£'bin,i +^bin,2) is liberated, either binary is 
destroyed, and its 1/4 goes to the kinetic energy of the center of mass of the surviving binary, and its 3/8 goes to each 
single star which is a component of the destroyed binary. Then, H{xi,X2) is expressed as 



H{xi,X2) 



< T 



2m|<I>r 



if {xi +X2) < 
16m|$ 



16m|$c 



4m|<i>c 
A^bb 
and substituting m|(f>c| = 

1, 



A-Bbb 



if (2:1 



kT 
fa;2)>16 



3 kT 

< (a;i+a;2) < 16 



to|$c 
kT 



(25) 



kT 



H{xi,X2) 



Xi +X2 

80 



(.Tl +X2) 



lOkT we finally obtain 
if (.Tl +2:2) < 53 ; 
if 53<(.Ti+a;2)<160; 

if (xi 4-2:2) > 160. 



40 



(26) 



Wc show the functions g'{x) and g{x)h{x) in figure 12, and the function G{xi,X2)H{xi,X2) in figure 13. These func- 
tions are required for calculating the dimensionless efficiency factors, expressed as equation (12) and (13). Additionally, 
we show the function g{x) in figure 12 and the functions G(2:i,a;2) in figure 13. In figure 12, the g{x)h{x) and g{x) 
curves overlap each other in 2: < 38. In figure 13, the G{xi,X2)H{xi,X2) curves and G'(xi,a;2) lines in Xi =X2, xi = 10x2, 
and xi ~ 1002:2 overlap each other in 2:2 < 0.5, 2:2 < 5, and 2:2 < 25, respectively. Note that the dimensionless heating 
rate of binary-binary encounters is not correct if 2:2 is soft, i.e. 2:2 < 3kT. This is because the cross section of the 
binary-binary interaction a is applicable only when the two binaries are hard. 
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3.3.5. Core size at the halt of core contraction 

We compare the core sizes of the clusters in our simulations with those derived from equation (10). When we use 
equation (10) to derive the ratio rc/r^ at the halt of core contraction, we approximate the distribution function of the 
binding energies /(x) and the unit of the binding energy kT as f{x) ~ d{x — i?bin,o/fc2^) and kT ~ A:To, respectively. 
The theoretical estimate of the core sizes is expressed as 

0.0196A fv^.rV , , ^ .. „ ,4 



rh logio(0.4Ar) 

/ 2/b.c(l-/b,c) 



(^) (^)<-"' 

[giEi,n,^o/kT)h{Eunfi/kT) - .g'(Sbi„.o/fcr)] 



+ l(j^^) G{Ei,infi/kT,Eun,o/kT)H{Ei,in,o/kT,Eunfi/kT)\. (27) 

We justify the approximation that /(a;)~(5(a; — £'bin.o/^'2^) a-nd kTc:ikTQ as follows. Figure 14 shows the distributions 
of the binding energies of binaries in the clusters of models lOkTo — 0.1, SOfcTo — 0.1, IOO/sTq — 0.1, and SOO/sTq — 0.1 at 
the time indicated in these panels, i.e. the time when the core contractions stop. All the distributions of the binding 
energies have steep peaks at the initial binding energies, i?bin,o- Even inside the core radii, all the distributions of the 
binding energies may do so. The distributions of the binding energies of binaries in the core at the halt of the core 
contraction can be also regarded as the same as the initial distributions, i.e. delta functions. 

The average kinetic energies of single stars within the half-mass radii and of binaries within the core radii are not so 
different from the initial time to the time at the halt of the core contraction, as seen in figure 9 and 10, respectively. 
We can regard IkT^lkTo. 

In figure 15, the big black dots show Tc/rh at the halts of core collapse of the clusters obtained in our simulations, 
as a function of the initial binding energy in the unit of kT, £'bin,o/^'2^- The numbers beside the dots show the mass 
fraction of the binaries in the core, /b,ci at the time of the halt, obtained from figure 8. The error bars indicate the 
amplitude of gravothermal oscillations. The dots are the geometric means of the maximum and minimum rc/rh in the 
gravothermal oscillations. Solid curves show equation (27) when /b,c = 0.04,0.1,0.4,1, which is shown by the numbers 
in italic format beside the curves. Here, we adopt Vg^/vh — \/2, 7 = 10, ant the free parameter A = 0.45. In the binding 
energy with less than E'bin ^ ^kin,avc.s ^ 2.5fcTo, the curves are not reliable, and should fall down to zero. Figure 15 
shows that the pairs of the ratios r^jr^^ and the binary fraction in the core /b,c at the halts of core contraction for 
models lOfcTg — 0.1, SOfcTo — 0.1, lOOfcTo — 0.1, and SOOfcTo — 0.1 in our simulations are included in those expressed by 
equation (27). 

The pairs of the ratios r^/ri^ and the binary fraction in the core /b,c at the halts of core contraction for models 
IfcTo — 0.1 and SfcTo — 0.1 arc not included in those expressed by equation (27), and the ratios rc/rh are much larger 
than equation (27) despite of their small binary fraction in the cores. This is because the core collapses of the clusters 
stop due to the energy that goes to the cluster from the three-body binaries composed of the single mass stars. The 
three-body binaries appear after the deep core collapse as discussed in section 3.2. 

We compare the pairs of the ratios rc/rh and the binary fraction in the core /b,c at the halts of core contraction 
predicted by equation (27) with those for model Double as the hard limit of the primordial binaries. Equation (27) 
should not include the simulation results of model Double. Equation (27) indicates that the ratio Tc/rh is zero, 
regardless of the mass fraction of the primordial binaries in the core, /b,c (see solid curves in figure 15). On the other 
hand, J'c/rh = 0.005 and /b,c ^ 1 at the halt of the core contraction in model Double. 

The reason for such disagreement is as follows. In the theoretical estimate which derives equation (27), only 
primordial binaries are considered as energy sources for the cluster. However, in model Double, another energy source 
appears. The energy source is binaries consisting of two double mass stars. The binaries are formed through the 
encounters of three single stars with double mass. The thick curve in model Double in figure 4 shows that such 
binaries generate energy. 

We expect that, in cluster models with harder primordial binaries than those we treat, similar energy sources to 
binaries consisting of two double mass stars in model Double appear. The energy sources arc hierarchical quadruple 
systems in which two binaries orbit around each other. The hierarchical quadruple systems may generate energy by 
shrinking orbits of the two binaries through interactions with the surrounding stars. 

We estimate the critical binding energy of the primordial binaries in which the hierarchical quadruple systems are 
formed, generate energy, and stop the contraction of the cluster core. We expect that the hierarchical quadruple 
systems are formed when rc/r-^ ~ 0.005. This is because the binaries consisting of two double mass stars are formed 
at such core size in model Double. When E'bin ^ lOOOfcToj the core contraction stops at rc/rh ^ 0.005 as seen in solid 
lines of /b.c = 0.4 and 1.0 in figure 15. It is not necessary to consider low mass fraction of the primordial binaries in 
the core, such as /b,c < 0.1, since the primordial binaries become more centrally-concentrated as they becomes harder 
(see figure 8). Therefore, the critical binding energy is ~ IOOO/cTq. 
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3.4- Initial mass fraction 

In this section, we discuss the dependence of the core evolution on the initial mass fraction of the primordial binaries, 
/b,o- Figure 16 shows the time evolution of the core radii, Tc, and half-mass radii, rh, of the clusters with /b,o = 0.03, 
0.1, and 0.3 primordial binaries, each of which has the binding energy E'bin.o = 3fcTo, 30fcro, and 300A:To. In models 
3kTo — 0.03, and 3fcro — 0.1, deep core collapse occurs, and in model 3kTo — 0.3, core collapse stops halfway. In all 
£'bin,o = 30A:To models, core collapse stops halfway. In all i^bin.o ~ 300fcTo models, deep core collapse occurs. 

Figure 17 shows the time evolution of the mass fraction of the binaries inside the core and half-mass radii, /b, 
of the clusters in the /b,o = 0.03, 0.1, and 0.3 models, each of which has the binding energy £'bin,o = 3fcTo, 30fcTo, 
and SOOfcTo. In all of them, mass segregation occurs initially. However, the mass fraction of the binaries in the core 
decrease halfway in models 3fcTo - 0.03, 3fcTo - 0.1, and 30fcro - 0.03. In models 30fcTo - 0.03, and 300fcTo - 0.03, the 
turning points correspond to the time when core collapse stops. In contrast, in models SkTo — 0.03, and SfcTo — 0.1, 
the mass fractions of the binaries in the core decrease long before the core collapse. 

In contrast to models 3fcTo — 0.03, and 3fcTo — 0.1, the core collapse of the cluster in model SfcTo — 0.3 stops halfway. 
Since the cluster has many primordial binaries, the energy from the primordial binaries to the cluster is large enough 
to stop the core collapse. 

We compare the ratios of the core radii to the half-mass radii in our simulation results with the theoretically 
estimated ratio in equation (10). The dots in figure 18 show the mass fraction of the binaries in the core at the halt 
of core collapse, /b,c, of the clusters whose primordial binaries have the initial binding energy, -Ebin.o- The numbers 
beside the dots show the ratios of the core radii to the half-mass radii at the halt of core collapse. When gravothermal 
oscillations occur, the geometric means are shown. The triangles, circles, and squares show the models /b.o = 0.03, 
0.1, and 0.3, respectively. 

The curves in figure 18 show the ratios of the core radii to the half-mass radii, Tc/ni, of clusters whose primordial 
binaries have the initial binding energy, i^bin.Oi and whose mass fraction of the binaries in the core is /b,c at the halt 
of the core collapse. The numbers in italic format beside the curves indicate the values of Tc/rh- If fi'bin.o ^ 2.5fcTo, 
the curves are not reliable. This is the same reason as in figure 15. 

In the clusters above the curve of rc/rh = 0.002, the core collapse stops halfway, and in the clusters below the 
curve of rc/rh = 0.002, the clusters experience deep core collapse. In the models -Bbin^o = SfcTo, model 3fcTo — 0.03 and 
SfcTo — 0.1 is below the curve, and model 3fcTo — 0.3 is above the curve. This is in good agreement with our simulation 
results. The ratios rc/rh in models IfcTg — 0.1, SfcTo — 0.03, and SfcTo — 0.1 disagree with theoretical curves, since the 
core collapse stops due to energy heating core generated by the three-body binaries. 

3.5. High-velocity escapers 

We investigate escapers of each cluster in the /b.o = 0.1 models. Figure 19 show the number of single escapers 
(iVosc,sin: solid lines) and binary escapers (iVcsc.bin^ dotted lines) in each logarithmic bin of velocities of the single 
escapers (usin) and the binary escapers (wbin) for all /b,o = 0.1 models. The velocities of these escapers are measured 
at the moment when they satisfy the conditions of the escapers, as shown in section 3.2. Note that in model Double, 
the single and binary escapers correspond to the escapers of the single and double mass stars, respectively. In the 
single escapers for all models except model Double, two peaks are present, although the higher peaks are small in 
models No-binary, IfcTg — 0.1, and SfcTg — 0.1. The population of the lower velocity escapers is driven by two-body 
relaxation, and that of the higher velocity escapers is ejected from the clusters through binary-single and binary-binary 
encounters. In the binary escapers, the escape velocities are similar to high-velocity population of the single escapers. 
They are also ejected from the clusters through binary-single and binary-binary encounters. 

Figure 20 shows the largest and the top 1 per cent (O.OlA'esc.sin-th largest), top 10 per cent (O.lA'csc.sin-th largest), 
and top 50 per cent (0.5A'osc,sin-th largest) velocities of single escapers from top to bottom in as a function of the 
initial binding energy of the primordial binaries, Ebin,o, in the models /b,o = 0.1. Figure 21 shows those of the binary 
escapers, although the top 1 per cent (O.OlA'cscbin-th largest)velocity is omitted because of the small number of the 
binary escapers. The dashed lines in both figures show circular velocities of the binaries as a function of the initial 
binding energy, £'bin,o- 

Consider a globular cluster whose virial radius is 10 pc, and whose mass is 10^ solar mass. Then, one velocity unit 
is 24:kms~^. When we apply our simulation results for the cluster, the highest velocity of the escapers is 500kms~^ , 
which is the single escaper in model 300fcro~0.1. In order to form hyper- velocity stars, which orbit in our Galaxy 
at speeds of 500- lOOOfcms^i (Hirsch et al. 2005; Brown et al. 2005; Edelmann et al. 2005; Heber et al. 2008), the 
globular clusters have to contain a large fraction of hard binaries with ^ 300A;To at the initial time. The presence of 
many hard binaries with '^ 300A:To at the initial time is possible, since SOO/cTq binaries are contact binary when the 
binary components are main-sequence stars. 

Figure 22 shows the number of binary escapers (A^oscbin) in each logarithmic bin of binding energies of the binary 
escapers (E'bin) for model No-binary and all /b,o ==0.1 models. Although the total numbers of the binary escapers arc 
different among model No-binary and all /b,o = 0.1 models, the binding energies of most binary escapers ranges from 
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lOOfcTo to lOOOfcTo. Figure 23 shows the largest, and the top 10 per cent (O.lA^csc.bin-th largest) and top 50 per cent 
(0.57Vcsc,bin-th largest) binding energies of the binary escapers from top to bottom as a function of the initial binding 
energies, i?bin,o for all /b,o = 0.1 models. 

Additionally, we list the triple escapers in table 3. The second column is the velocity of the center of mass of the 
triple escapers, Wtri- The third and fourth columns are, respectively, the binding energy (i^bin.in and Sbin.out) of the 
inner and outer binaries in unit of fcTp. 

4. Summary 

We study systematically the dependence of cluster evolution on the binding energy of primordial binaries. By means 
of GORILLA, we simulate the core evolution of the clusters, each of which contains primordial binaries with equal 
binding energy. 

When the initial mass fraction of the primordial binaries is fixed to 0.1, we find that the dynamical evolutions of 
the clusters are divided into three ranges according to hardness of the primordial binaries as follows. 

1. In soft range (< SfcTg), the clusters experience core collapse in similar way to those without primordial binaries. 
The ratios of core radii to half-mass radii at the halt of the core collapse are about 0.006. The primordial binaries 
do not heat the clusters. This is because the primordial binaries arc destroyed through encounters with single 
stars, and do not generate energy. 

2. In intermediate hard range (lOfcTo — IOO/cTq), the core collapses in the clusters halt halfway. The ratios of core 
radii to half-mass radii at the halt of the core collapse are 0.05 — 0.1. The primordial binaries release energy, 
and the energy heats the clusters. 

3. In super hard range (> 300A:To), the clusters experience core collapse, and the ratios of core radii to half-mass 
radii at the halt of the core collapse is about 0.02. The primordial binaries do not so much heat the clusters. 
Although the primordial binaries release energy through encounters, the energy is so large that binaries and 
single stars involved with the encounters are ejected from the clusters. 

The dependences of the boundaries between the soft and intermediate hard ranges and between the intermediate 
and super hard ranges on the initial mass fraction of the primordial binaries are as follows. 

1. The boundary between the soft and intermediate hard ranges depends on the initial mass fraction of the pri- 
mordial binaries. When the mass fraction of the primordial binaries is 0.3, the core contraction in the cluster 
with SkTo primordial binaries halts at large ratio of core radius to half-mass radius ^ 0.07, and the intermediate 
hard range includes SkTo. 

2. The boundary between the intermediate and super hard ranges is not changed, when the initial mass fraction of 
the primordial binaries ranges from 0.03 to 0.3. 

We compared the pairs of the ratios of core radii to half-mass radii and the core mass fraction of the binaries at 
the halt of the core contraction in our simulations with those of theoretical estimates. We found a good agreement 
between iV-body simulations and the theoretical values. 
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Table 1. Initial models. 



Model name 



E, 



bin,0 



/b,0 



N, 



b.O 



/d.O 



N, 



d.O 



IfcTo- 


0.1 


IfcTo 


0.1 


819 








SfcTo- 


0.1 


SfcTo 


0.1 


819 








10A:To - 


-0.1 


lOfcTo 


0.1 


819 








30/cTo - 


-0.1 


SOfcTo 


0.1 


819 








lOOfcTo 


-0.1 


100/cTo 


0.1 


819 








SOOfcTo 


-0.1 


300/sTo 


0.1 


819 









SfcTo - 0.03 
SOfcTo - 0.03 
300A:ro - 0.03 


3fcro 

30fcro 

300A:ro 


0.03 
0.03 
0.03 


246 
246 
246 












3fcro - 0.3 
30A:To - 0.3 


ikTo 
30kTo 


0.3 
0.3 


2458 
2458 










No-binary 
Double 













0.1 



819 



Table 2. Accuracy, and apoccntric and pericentric parameters. 



Model name 


V 


Vs 


a 


/3 


lfcTo-0.1 


0.01 


0.0025 


5 


10 


3fcro-o.i 


0.01 


0.0025 


5 


10 


lOkTo-0.1 


0.01 


0.0025 


5 


10 


30/cTo-O.l 


0.01 


0.0025 


5 


10 


lOO/cTo-O.l 


0.01 


0.0025 


5 


50 


300A:To-0.1 


0.01 


0.0025 


8 


50 


3kTo - 0.03 


0.01 


0.0025 


5 


10 


30fcro - 0.03 


0.01 


0.0025 


5 


10 


300fcTo - 0.03 


0.01 


0.0025 


5 


50 


3fcTo - 0.3 


0.01 


0.0025 


5 


10 


30A:To - 0.3 


0.01 


0.0025 


5 


50 


No-binary 


0.01 


0.0025 


5 


10 


Double 


0.01 


0.0025 


5 


10 



Table 3. The list of triple escapers. 



Model name iJtri [standard units] -Bbin.ini^To] -E-bin.out [^To] 



30fcro-0.1 4.7 



8.5 X 10^ 
2.8 X 10^ 
5.7 X 10^ 
4.7 X 10^ 
6.2 X 10^ 
9.4 X 10^ 



3.0 X 10"^ 
9.2 X 10" 

7.4 X 10° 
2.6 X 10° 

1.5 X 10° 
3.7x10^ 



lOO/cTo-O.l 



0.44 
1.2 
0.39 
1.7 



300fcro-0.1 2.0 
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Fig. 1. Illustration of binaries in isolation with conditions (A) (upper panel), and conditions (B) (lower panel). 
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Fig. 2. Energy errors as a function of simulation time. 
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Fig. 3. Time evolution of the core radii, Tc, and half-mass radii, r'h, of /b,o = 0.1 cluster models, and models No-binary and Double. 
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30202020 
t (x 1000) 



Fig. 4. The increase of the total binding energy of the binaries, ABbin.totC*) (thick curves), which corresponds to energy generated 
by the binaries, and the total kinetic energy of escapers, -Eesc,tot(t), from the clusters (thin curves). The arrows indicate the times 
when the core collapse stops. 
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Fig. 5. Time evolution of the number of binaries, A'^^i in the /^.o = 0.1 models and models No-binary and Double. For model 
Double, the number of the double mass stars is also plotted. The thick curves indicate the numbers of binaries (or double mass 
stars) within the clusters, and the thin curves indicate the total numbers of binaries, (or double mass stars) including escapers. The 
arrows indicate the times when the core collapse stops. 
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Fig. 6. The number of binaries, A'^^i iii each logarithmic bin of the binding energy, -Ebin 
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Fig. 7. Binding energies as a distance from cluster center of each binary at the time indicated in the panels. The dashed lines show 
the half-mass radii and twice the core radii at the time. 
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Fig. 8. Time evolution of the mass fraction, /b, of the binaries inside the core radii (the upper curves) and half- mass radii (the 
lower curves) in the /b o = 0-1 models and model No-binary. For model Double, the mass fraction, /j, of the double mass stars are 
shown. The arrows show the time when the core collapse stops. 
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Fig. 9. Time evolution of the mean kinetic energy of tlie single stars inside the core radii (solid curves) and the half-mass radii 
(dashed curves) in the /b o = 0-1 models and models No- binary and Double. The arrows show the time when the core collapse stops. 
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Fig. 10. Time evolution of the mean kinetic energy of the binaries, -Ekin.avc bi inside the core radii in the /^.o = 0.1 models. The 
arrows show the time when the core collapse stops. 
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Fig. 11. Interpretation of core evolution. 
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Fig. 12. Dimensionless heating rate of binary-single encounters in which the binaries survives (g{x)h{x)) and are destroyed {g'{x)) 
as a function of dimensionless binding energy x (sohd curve) . The dashed line shows the dimensionless hardening rates of the binary 
with the binding energy a: in a sea of single stars, i.e. g{x). 
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Fig. 13. Dimensionless heating rate of binary-binary encounters, G{xi,X2)H{xi,X2), as a function of X2, where the binaries have 
dimensionless binding energy xi and X2, and xi = X2, x\ = 10x2, and x\ = lOOxi (solid curves). The dashed lines show the 
dimensionless hardening rates of the binary with the binding energy xi in a sea of binaries with the binding energy X2, i.e. 
G{xi,X2). 
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Fig. 14. Distributions of the binding energies of binaries in the whole clusters at the time indicated in each panel, which is the time 
when the core contractions stop. 
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Fig. 15. Ratio of core radii to half-mass radii at the halts of core collapse of the clusters whose primordial binaries have equal 
binding energy, -Bbim a-nd whose cores contain the mass fraction of the primordial binaries in the core, /]-, (.• The dots show the 
ratio of the core radii to the half-mass radii at the halt of core collapse in our simulation. The numbers beside the dots are /]-, ^ at 
that time. The four curves draw equation (10) when /b c = 0.04, 0.1, 0.4, and 1.0. The values of /b ^ are beside the curves in italic 
format. The dashed lines, rc/rj, = 0.002, and 0.005, show the ratio of core radii to half-mass radii at the halts of core collapse in 
model No-binary and Double, respectively 
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Fig. 16. Time evolution of the core, re, and half-mass radii, rn, of the clusters with /b,o = 0.03, 0.1, and 0.3 primordial binaries, 
each of which has the binding energy i?bin = 3fcTo, lOfcTo, and 300fcTo. The way of calculation of the core radii is the same as that 
in figure 3. 
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Fig. 17. Time evolution of the mass fraction of the binaries inside the core and half-mass radii of the clusters, /b, with /b,o = 0.03, 
0.1, and 0.3 primordial binaries, each of which has the binding energy -Ebin,o = 3fcTo, lOfcTo, and SOOfcTg. 
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Fig. 18. Contours of the ratio of core radii to half-mass radii at the halt of core contraction in clusters whose mass fraction of 
primordial binaries in the core is /b,ci a-nd distribution function of the primordial binaries is 5{x — i?bin)i which are obtained from 
equation (27). The numbers in italic show the ratio of the core radii to the half-mass radii. The dashed curves are possibly not 
correct, since the heating rate through binary-binary interactions in hard binaries is extended to the soft range. The dots show 
simulation results. The vertical axis shows /t, ^ at the halt of core contraction, the horizontal axis shows Ebin a-t the initial time, 
and the numbers beside the dots is the ratio of core radii to the half-mass radii at the halt of the core contraction. The shapes of 
the dots show the initial mass fraction of the primordial binaries in the clusters, /b q. The triangles, circles, and squares show the 
models fho — 0.03, 0.1, and 0.3, respectively. 
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Fig. 19. The distribution of the velocities of single and binary cscapers (solid and dotted lines, respectively) in the models /j, g = 0.1. 
The velocity is in A^-body standard units. 
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Fig. 20. The largest, and the 0.01A^csc,sin-th, 0.1Afesc,sin-th, and 0.5A^csc,sin-th largest veloeities of single escapers from top to 
bottom as a function of initial binding energy of primordial binaries in the /j, o = 0.1 clusters. The dashed line shows circular 
velocity of a binary as a function of its binding energy. 
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Fig. 21. The largest, and the O.lA^csc bin-th and 0.5A^csc.bin-th largest velocities of binary escapers from top to bottom as a function 
of initial binding energy of primordial binaries in the /^^o = 0-1 clusters. The dashed line shows circular velocity of a binary as a 
function of its binding energy. 
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Fig. 22. The distribution of the binding energy of binary escapers in the models /b.o = 0.1. The binding energy is in the unit of 

kTo. 
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Fig. 23. The largest, and the O.lA^osc bin-th and O.SA^eacbin-th largest binding energies of binary escapers from top to bottom as a 
function of initial binding energy of primordial binaries in the /b,o = 0.1 clusters. 



